The generator coordinate method in time-dependent density-functional theory: 

memory made simple 



E. Orestes, 1 - 2 - 3 K. Capelle, 3 A. B. F. da Silva, 1 and C. A. Ullrich 2 

1 Departamento de Quimica e Fisica Molecular, Instituto de Quimica de Sao Carlos, 
Universidade de Sao Paulo, Caixa Postal 780, Sao Carlos, Sao Paulo 13560-970, Brazil 
2 Department of Physics and Astronomy, University of Missouri, Columbia, Missouri 65211, USA 
Departamento de Fisica e Informdtica, Instituto de Fisica de Sao Carlos, 
Universidade de Sao Paulo, Caixa Postal 369, Sao Carlos, Sao Paulo 13560-970, Brazil 

(Dated: February 1, 2008) 

The generator coordinate (GC) method is a variational approach to the quantum many-body prob- 
lem in which interacting many-body wave functions are constructed as superpositions of (generally 
nonorthogonal) eigenstates of auxiliary Hamiltonians containing a deformation parameter. This 
paper presents a time-dependent extension of the GC method as a new approach to improve exist- 
ing approximations of the exchange-correlation (XC) potential in time-dependent density-functional 
theory (TDDFT). The time-dependent GC method is shown to be a conceptually and computation- 
ally simple tool to build memory effects into any existing adiabatic XC potential. As an illustration, 
the method is applied to driven parametric oscillations of two interacting electrons in a harmonic 
potential (Hooke's atom). It is demonstrated that a proper choice of time-dependent generator 
coordinates in conjunction with the adiabatic local-density approximation reproduces the exact lin- 
ear and nonlinear two-electron dynamics quite accurately, including features associated with double 
excitations that cannot be captured by TDDFT in the adiabatic approximation. 



I. INTRODUCTION 

The development and applications of time-dependent 
density-functional theory (TDDFT) [J Q have enjoyed 
exponential growth in the last few years. The most suc- 
cessful applications have so far been in the linear or per- 
turbative regime, to calculate excitation energies and op- 
tical spectra of complex molecular systems [3| . Other ar- 
eas of interest for TDDFT are the nonlinear, nonpertur- 
bative electron dynamics of atoms and molecules driven 
by intense laser fields, and quantum transport through 
molecular junctions (see Ref. Q for a comprehensive re- 
view of formalism and applications of TDDFT) . 

The central quantity in TDDFT, the time-dependent 
exchange-correlation (XC) potential v xc (r, t), has to be 
approximated in practice, and the success (or the fail- 
ure) of a TDDFT calculation depends critically on the 
quality of the chosen XC functional. The most widely 
used approximation for the XC potential in TDDFT is 
the adiabatic approximation: 



adiabatic / 



(r,t) 
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where some given XC potential t£* atIC [n] from static DFT 
is evaluated at the time-dependent density n(r,t). This 
idea seems certainly reasonable for systems that are 
slowly varying in time. The most common example of 
this approach is the adiabatic local-density approxima- 
tion (ALDA) where, in addition to a slow change over 
time one also assumes slow spatial variations. 

The adiabatic approximation means that all functional 
dependence of v xc (r, t) on prior time-dependent densities 
n(r', t'), t' < t, is ignored. Neglecting the retardation im- 
plies frequency-independent and real XC kernels in lin- 
ear response [a, 0] . This approach is known to work well 



for excitation processes in many-body systems that have 
a direct counterpart in the Kohn-Sham system, such as 
atomic and molecular single-particle excitations. On the 
other hand, for more complicated processes such as dou- 
ble or charge-transfer excitations the ALDA can fail dra- 
matically @, H, [1] . We mention that the well-known fail- 
ure of the ALDA to reproduce optical absorption spectra 
in extended systems [l0| has less to do with the neglect 
of retardation, but is due to the wrong spatial long-range 
behavior of the linear-response XC kernel. 

Another well-known class of problems where the adia- 
batic approximation fails are correlated double photoion- 
ization processes in atoms and molecules (for a review, 
see [ll|). The most notorious case is the direct double 
ionization of helium [l2j], which has been a persistently 
nagging problem for TDDFT in any of the standard XC 
approximations [ll, E3, EH • 

Systematic studies of plasmon-like charge-density os- 
cillations in quantum wells and quantum rings [ToL E3, 
El, EU have uncovered further examples where it is nec- 
essary to go beyond the ALDA. Memory effects have been 
shown to be crucial in these examples, in particular when 
multiple excitations are involved. 

A fully-fledged treatment of memory and retardation 
in TDDFT presents many technical difficulties, although 
quite a bit of progress has been made recently. The most 
fruitful approaches so far are based on ideas of hydrody- 
namics and elasticity theory, where the memory resides 
in the dependence of a fluid element's motion and de- 
formation on its initial position, and where retardation 
effects cause the electron liquid to be subject to viscoelas- 
tic stresses [U [H J2I [2I HU . Applications of these 
ideas to polymers [26l \27\ and quantum wells [28| have 
met with some success, but there are also fundamental 
problems. A particularly disquieting situation arises in fi- 
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nite systems, where the hydrodynamic TDDFT approach 
results in excitation energies with unphysical imaginary 
parts [H, H • 

The purpose of this paper is to present a new, simple 
and intuitive approach to retardation effects in TDDFT. 
The approach is based on an extension the so-called gen- 
erator coordinate (GC) method, which was originally in- 
troduced in nuclear physics by Wheeler and coworkers 
[30l l3lj . The GC method expresses the many-body wave- 
function in terms of a superposition of states arising from 
deformed auxiliary Hamiltonians, using a variational op- 
timization of the expansion coefficients. In the present 
paper, we derive and apply a time-dependent general- 
ization of the GC method, combining it with TDDFT. 
The result is an approximation to the time-dependent 
many-body wave function, built from Slater determinants 
which come from deformed time-dependent Kohn-Sham 
Hamiltonians. We will show that by introducing a class 
of "deformations" involving the previous history of the 
system, a cheap and simple way to simulate retardation 
effects is obtained. 

The paper is organized as follows. Section [TT] first re- 
views the basic GC formalism, and then introduces the 
time-dependent GC method in the context of TDDFT. 
Our model system, Hooke's atom, is briefly discussed in 
Section HTT1 In Section ITVl we present applications of the 
GC-TDDFT approach to a simple model system (para- 
metric oscillations of two electrons in a harmonic poten- 
tial, also known as Hooke's atom). Section [V] gives a 
summary and conclusions. Some technical details are 
given in the appendices. We use atomic (Hartree) units 
everywhere, with e = m = K = 1. 



II. FORMAL FRAMEWORK 

A. The GC method 

In the GC method, approximations for the full many- 
body wave function ^(xi, ...,xjv) of a system of N elec- 
trons, where x = (r, s) denotes spatial coordinate and 
spin, are constructed as follows [3l| : 



*(xi 



,Xjv) 



daf(a)$(a,Tti, ...,xjv) 



(2) 



Here, $(a,Xi, ...,xjv) represents a set of (usually nonin- 
teracting) seed functions which depend on a deformation 
parameter a, the generator coordinate. This dependence 
can be explicit (if a is a parameter in an analytic expres- 
sion for or implicit (e.g., if the <£> are eigenstates of a 
Hamiltonian which contains the parameter a). 

The weight function f(ct) is determined variationally. 
Variation with respect to /* (a) of the normalized many- 
body expectation value where H is the 
Hamiltonian of the fully interacting system, leads to the 
so-called Griffin-Hill- Wheeler (GHW) equation: 



where K(a,a') = (&{a)\H |*(a')> and S(a,a') = 
($>(a)\$>(a')) are the Hamiltonian and overlap matrix el- 
ements of the given seed functions. Solution of the GHW 
equation © yields a set of energy eigenvalues E and as- 
sociated cigenfunctions f(a) which are used in Eq. ([2]) to 
determine the corresponding many-body wave functions 
and all observables following therefrom. 

In the original work of Wheeler et al. [3(1 HH , VP was a 
nuclear many-body wave function, and the generator co- 
ordinate was a parameter describing the deformed shape 
of the potential due to the collective behavior of the nu- 
cleons. The seed functions were obtained from a simpli- 
fied Schrodinger-like equation featuring a deformed nu- 
clear potential. 

Over the years, the GC method has found widespread 
use in a large variety of problems in nuclear and electronic 
many-body physics |32l . l33l l34j . In quantum chemistry, 
the GHW variational equations were used to develop 
high-quality Hartree-Fock and Dirac-Fock single-particle 
basis functions [H, HI) H3) HI]- Expressions similar to 
the GC ansatz ^ have recently also been used by Alon, 
Streltsov and Cederbaum [39j in order to build symme- 
tries into the many-body wave function \E' that are not 
present in the seed wave functions $; and by Pan, Sahni 
and Massa [I(| with the aim to allow more general varia- 
tions of the wave function during variational calculations. 



B. Static GC-DFT 

Recently, a connection between the GC method and 
static DFT was established by Capelle 



41] 



The idea is 

to express the electronic many-body wave function @ 
as a weighted superposition of Kohn-Sham (KS) Slater 
determinants $KSi constructed with KS orbitals ip": 



*(xi, ...,xjv) = / da/(a)$ K s(o;,xi, ...,xjv). 



(4) 



The generator coordinate is introduced, like in the origi- 
nal proposal of Wheeler et al. [13, [H| , as a deformation 
parameter of the potential, in this case of the XC poten- 
tial. The KS orbitals (p? , j = 1, . . . ,N, following from 
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+ V ext (r) +v H (r) +<.(r) 



p?(r) = e>?(r) , 



da [K(a, a') — ES(a, a')) f(a') = 0, 



(3) 



(5) 

where w ex t an, 4 v u are the external and Hartree poten- 
tials, used to construct the Hamiltonian and overlap ma- 
trix elements K(a, a') and S(a, a') in the GHW equation 
Q. In practice, the GC-DFT wave function ^ is con- 
structed using a finite set of deformation parameters a 
(usually less than 10), forming a mesh of points which is 
used to discretize the GHW integral equation ([3]). Here 
and in the following, we only consider systems that are 
nonmagnetic everywhere and thus ignore the spin index. 

A first illustrative numerical calculation of this new 
approach, showing its viability and potential, was made 
in Ref. [4l| for the Helium isoelectronic series. In this 



3 



case, a single set of deformation parameters was tested 
together with a simple functional, the X-only LDA, to 
obtain ground-state energies and the weight functions. 
Through Eq. ^ an approximation to the many-body 
wave function was established and the expectation val- 
ues of the operator r n , with n = —2,-1,0,1,2, were 
calculated. 

The results of Ref. [4l| show that the ground-state 
energies obtained are considerably better than those cal- 
culated from both the LDA generator functionals them- 
selves and more sophisticated gradient-corrected XC 
functionals. The many-body wave functions were ob- 
tained with almost no additional numerical effort, yield- 
ing expectation values that are close to those given by 
other methods. Thus, the GC method is a cheap and 
conceptually simple way to improve existing XC func- 
tionals. The method is, of course, not restricted to the 
helium isoelectronic series. Moreover, it can also be used 
to extract excitation energies. Work along these lines will 
be reported elsewhere [42| . 



C. Time-dependent GC-DFT 

Let us now extend the GC-DFT approach into the time 
domain. Our starting point is a generalization of Eq. (|2j): 
we write the time-dependent many-body wave function 'J 
as a linear superposition of time-dependent KS (TDKS) 
Slater determinants <&ks depending on the generator co- 
ordinate a: 



*(xi, ...,x N ,t) = J rfa/(a,t)$Ks(a,xi, . 



,x N ,t). (6) 



The weight function /(a, t) is taken to be time-dependent 
and complex. We consider situations where the system 
is in a stationary state for t < t$. The time evolution 
of \I/ for t > to thus represents an initial value problem, 
where ^P, $ks an d / at the initial time to are obtained 
from the static GC-DFT method described above. 

Like in the static case, the KS Slater determinants are 
built from single-particle orbitals 0"(r, t), j = 1, . . . ,N, 
where 



V 2 d 

— + Wext (r, t) + v H (r, t) + w° (r, t) - i— 



0?(r,t)=O. 
(7) 

The initial conditions are </>"(r, to) = tp" (r) for each a, 
see equation §5$. 

The time-dependent weight functions f(a,t) are deter- 
mined through a stationary-action principle: 



5f*(a,t) 



dt' ( (*') 



,d_ 



Hit') 



0, 



(8) 

where t\ is an arbitrary upper limit of the time propa- 
gation interval. H(t) is the time-dependent many-body 
Hamiltonian featuring the external one-body potential 
■u ex t(r, t). 



The solution of Eq. ([8]) leads to the time-dependent 
GHW (TDGHW) equation: 



da' 



A(a, a' , t) + S(a, a' , t)i 



where 
A(a, a' , t) = 



$ks(">*) 



dt 



H{t) 



/(o',t) = 0, (9) 



*Ks(a',t)) (10) 



are the KS action matrix elements, and 

S(a,a',t) = (<S>* KS (a,t)\<I> KS (a',t)) 



(11) 



are the time-dependent KS overlap matrix elements. In 
Appendix [X] we show that the TDGHW equation sat- 
isfies two key properties: it reduces to the static GHW 
equation ([3]) in the limit where v ex t is time- independent, 
and it conserves the norm of the total wavefunction ^ . 

In practice, one first chooses a mesh of a- values and ob- 
tains the initial KS orbitals (r, to) and the weight func- 
tions f(a,to) from a static GC-DFT calculation. One 
then propagates the TDKS equation ([7]) for each KS or- 
bital and each value of a, from to to t\. Then, the ma- 
trix elements ([10} and (|TTJ) are formed, and the TDGHW 
equation can be solved by time propagation, which yields 
the time-dependent weight functions f(a,t). The final 
step is to calculate the observables of the system under 
study. 

In Appendix |5] we discuss our numerical approach to 
solve the discretized TDGHW equation using the Crank- 
Nicholson algorithm, which guarantees a unitary time 
evolution of 'l'(t). 



III. HOOKE'S ATOM: MODEL AND 
OBSERVABLES 

To illustrate the GC-TDDFT approach, we consider a 
system of two interacting electrons confined by a spher- 
ically symmetric harmonic potential kr 2 /2, known as 
Hooke's atom. The great advantage of this system is 
that the associated two-body Schrodinger equation sep- 
arates into center-of-mass and relative coordinates, and 
can therefore be solved in a straightforward manner, even 
in the time-dependent case [13, S3, HE, S^, S3, IH, Ha, H3] - 
This allows us to compare our approximate TDDFT and 
GC-TDDFT results with exact numerical solutions of the 
full Schrodinger equation, at a relatively small computa- 
tional cost. 

In the following, we shall consider Hooke's atom driven 
by a time-dependent spring constant 



k(t) = [1 + As(t) sin(w<)]fc 



(12) 



where k is the spring constant of the initial ground state, 
A and u> are the amplitude and frequency of the oscilla- 
tions, and s(t) defines a "pulse shape". 
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Hooke's atom resembles the helium atom in many re- 
spects, but there is one important difference: due to the 
harmonic confining potential, Hooke's atom has no con- 
tinuum of unbound states and hence cannot be ionized. 
Therefore it cannot be used to address the issue of direct 
double ionization [T3, EH EH EH . Nevertheless, Hooke's 
atom is a useful model system for studying correlated 
electron dynamics, including double excitations. 

It is known that TDDFT provides the time-dependent 
density n(r, t) in principle exactly, but runs into difficul- 
ties whenever one is interested in observables that are 
not easily expressed as explicit functionals of the density. 
Thus, a quantity such as the radial expectation value 



(r(t)> 



d 3 r n(r, t) 



(13) 



can be easily calculated, since it involves only the total 
density. By contrast consider the following quantities for 
a two-electron system [U, HH : 

P (0} (t) = \ f d 3 n f ^^(n,!-^)! 2 ^) 

1 Jri<R Jr 2 <R 

P (2 \t) = \( d 3 n f rf 3 r 2 |*(n,r 2 ,t)| 2 (15) 

Z Jr x >R Jr 2 >R 

pW(i) = 1-Pt°)(i)-P( 2 )(t), (16) 

where P(°) and P( 2 > are the probabilities of finding both 
electrons inside or outside of an imaginary spherical box 
of radius R surrounding the system, and P^ is the prob- 
ability of finding one electron inside and the other one 
outside the box. For a two-electron real atom, these 
quantities are ionization probabilities, provided R is large 
enough. 

By the Runge-Gross theorem Q, P<-°\ pW and P^ 
have the formal property of being functionals of the time- 
dependent density. The explicit form of these functionals, 
however, is unknown, and it is therefore very difficult to 
extract 

p(o) i p(i) and p(2) only from n ( r ,i). It is com- 
mon practice (although not formally justified) to replace 
the full interacting wave function ^ in Eqs. (jUJ)— (|T5J) 
by the KS wave function <J>ks [13 • However, it comes 
as no surprise that the so-defined KS probabilities fail to 
capture the essence of the correlated double ionization 
processes in helium [l3|. A method that is wavefunction 
based, such as GC-DFT, has a much better chance for 
success. 
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FIG. 1: (r(t)) for a parametrically driven Hooke's atom, with 
spring constant k(t) given by Eq. I|12jl. Top and bottom 
panels: A — 0.1 and 0.25; left and right panels: uo = 1.9 and 
2.0 (above and below resonance). The "pulse shape" of k(t) 
is trapezoidal, with a duration of 9 cycles. Solid lines: exact 
calculation. Dashed lines: ALDA. 
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FIG. 2: Probability P (0) [Eq. fJH}] of finding both electrons 
inside a spherical box of radius R = 2.7 a.u., for the system of 
Fig. [1] Solid lines: exact calculation. Dashed lines: ALDA. 



at constant amplitude A, and a two-cycle linear switch- 
off. Driven by the time-dependent spring constant, the 
system carries out breathing-mode oscillations which pre- 
serve the initial spherical symmetry. We shall consider 
driving frequencies u> in the vicinity of the lowest reso- 
nance at ujq — 1.95 a.u., for fco = 1 a.u. 



IV. RESULTS AND DISCUSSION 



ALDA 



We now present results for parametrically driven oscil- 
lations of Hooke's atom, comparing exact solutions of the 
time-dependent two-particle Schrodinger equations with 
various levels of TDDFT and time-dependent GC-DFT. 

In the following, we use time-dependent spring con- 
stants k(t) with a trapezoidal shape for s(t), with a 
switch-on by a linear two-cycle ramp, followed by 5 cycles 



We first compare (x)-only ALDA with exact results 
for the time-dependent Hooke's atom. Fig. [1] shows 
(r(t)), and Figs. HHSpresent p(°\ pW and P {2 \ respec- 
tively, for four different combinations of amplitudes and 
frequencies of the driving force constant (A — 0.1,0.25, 
and us — 1.9,2.0). The spherical box used to compute 
the probabilities has radius R — 2.7 a.u. 
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FIG. 3: Probability P (1) [Eq. JTBJ] of finding one electrons 
inside a spherical box of radius R = 2.7 a.u., and the other 
one outside, for the system of Fig. [T] Solid lines: exact 
calculation. Dashed lines: ALDA. 
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FIG. 4: Probability P (2) [Eq. JEJ] of finding both electrons 
outside a spherical box of radius R = 2.7 a.u., for the system 
of Fig. [1] Solid lines: exact calculation. Dashed lines: ALDA. 



The ALDA completely fails to reproduce the behav- 
ior of the driven Hooke's atom. This can be seen very 
clearly from Fig. Q] in ALDA, the radial expectation 
value oscillates with a constant amplitude after the trape- 
zoidal pulse is over. By contrast, the exact (r(t)) exhibits 
a beating pattern, where the amplitude of the oscilla- 
tions is strongly modulated. Likewise, the probabilities 
p(o) j pi 1 ) and P^) show a drastic amplitude modulation, 
which is absent in ALDA. 

The physical origin of the beating patterns in the 
charge-density oscillations is a superposition of the dy- 
namics associated with single and multiple excitations. 
For Hooke's atom, the beating pattern is particularly 
pronounced since, due to the harmonic confining poten- 
tial, the electronic level separations are nearly equidis- 
tant, which makes the coupling to double and higher 
excitations very efficient. A similar though less pro- 
nounced modulation effect of charge-density oscillations 
was recently analyzed for two electrons confined on a 
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FIG. 5: Same as Fig. [T] Solid lines: exact calculation. 
Dashed lines: adiabatic GC-TDDFT. 
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FIG. 6: Same as Fig. [2] Solid lines: exact calculation. 
Dashed lines: adiabatic GC-TDDFT. 



two-dimensional quantum strip (l9( . 

It is a well-known fact that the ALDA fails to capture 
any type of excitations which does not have a counter- 
part in the Kohn-Sham single-particle spectrum, such as 
double excitations [8, 9]. Nevertheless, the sheer magni- 
tude of the error of ALDA exhibited in Figs. [THU comes 
somewhat as a surprise. 



B. Adiabatic GC-TDDFT 

In a first attempt to improve the ALDA description 
with the GC-TDDFT approach, we consider the follow- 
ing deformed xc potential: 



i£ c (r,t) = (a - l)v a (r,t) + a<4 UA [n(r,*)] 



(17) 



For a = 1, this reduces to the ALDA, and for a = 
the xc potential vanishes and the first term cancels the 
Hartree potential in the TDKS equation, leaving only 
the bare confining potential. This describes the situation 
where one electron is far away from the center, and the 
other electron only sees the bare potential. Due to the 
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FIG. 7: Same as Fig. [3] Solid lines: exact calculation. 
Dashed lines: adiabatic GC-TDDFT. 




FIG. 9: Same as Fig. \T\ Solid lines: exact calculation. 
Dashed lines: nonadiabatic GC-TDDFT. 
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FIG. 8: Same as Fig. [4] Solid lines: exact calculation. 
Dashed lines: adiabatic GC-TDDFT. 
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FIG. 10: Same as Fig. [2] Solid lines: exact calculation. 
Dashed lines: nonadiabatic GC-TDDFT. 



C. Nonadiabatic GC-TDDFT: memory 



self-interaction error, the ALDA does not capture such 
a scenario correctly, whereas the deformed XC potential 
(fTT)) . for a — 0, cancels this error. 

We have chosen the set of deformation parameters 
a = {0, 1, 2, 3, 4} and the x-only ALDA in Eq. ([T7J). The 
results of this adiabatic GC-TDDFT approach for the 
charge-density oscillations in Hooke's atom are shown in 
Figs. [SHU f° r the same parameters as in Figs. HHU We 
find that the improvement of the adiabatic GC-TDDFT 
over the pure ALDA is marginal. There is a hint of an 
amplitude modulation in all cases, but it is by far insuf- 
ficient compared to the exact data. 

A possible reason for the failure of this method is the 
peculiar nature of Hooke's atom which prevents ioniza- 
tion. It seems plausible that in a real atom, where one 
or both electrons can ionize, the approximate treatment 
of the self-interaction error implied in Eq. (|17p might be 
more beneficial. 



As discussed in the introduction, the adiabatic ap- 
proximation omits memory and retardation effects, which 
have been shown to be crucial for the treatment of dou- 
ble and multiple excitations @. Therefore, we will now 
choose a deformation in the GC-TDDFT approach which 
introduces retardation in a very simple manner. Let 

t&(r,t) - (a- lK(r,i) + av^ A [n(r,t r )} , (18) 

where the retarded time, t r , is given by 

aT 

U=t- — . (19) 

Here, T = 2tt/u> denotes one cycle of the time-dependent 
spring constant used to trigger the charge-density oscil- 
lations. Equations (|18[) and (flU)) imply that the LDA xc 
functional at time t is evaluated using the time-dependent 
density at a previous time t r , where a determines how 
far we go back into the past. Choosing again the set of 
deformation parameters a = {0,1,2,3,4}, we go back 
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A=0.10andcu=1.9 



mmw 1 



A=0.1 and m=2.0 



A=0.1 and <d=2.0 



■■IIP 1 




A=0.1 and a)=2.0 




20 40 60 80 100 20 -10 60 80 100 

time (fs) time (fs) 

FIG. 11: Same as Fig. [3] Solid lines: exact calculation. 
Dashed lines: nonadiabatic GC-TDDFT. 




FIG. 12: Same as Fig. [4] Solid lines: exact calculation. 
Dashed lines: nonadiabatic GC-TDDFT. 



at most half a cycle, which corresponds to a phase lag 
of 180° of the deformed xc potential with respect to the 
charge-density oscillations. 

Results are shown in Figs. I9TTT21 for the same param- 
eters as in Figs. [5H5] and [THU We obtain an excellent 
agreement between the exact and the nonadiabatic GC- 
TDDFT results for (r(i)), where the strong amplitude 
modulations of the charge-density oscillations are very 
well reproduced for all four parameter combinations in 
Fig. The improvement for the probabilities ; Figs. 
ITOt4"T2l is also very substantial, especially for the larger 
amplitude A = 0.25. 



D. Weight Functions 

In order to analyze the contribution of each TDKS 
determinant in the approximation of the time-dependent 
many-body wave function, ^(ri, rjy, t), we plot the 
square modulus of the weight functions \f(a,t)\ 2 in Fig. 
HUI for the nonadiabatic GC-TDDFT calculation with 
A = 0.1 and u> = 1.9. We will focus in particular on the 




20 10 

time (fs) 



FIG. 13: Weight functions \f(a, t)\ 2 for the nonadiabatic GC- 
TDDFT calculation of the driven Hooke's atom, with A — 0.1 
and lo = 1.9. 



two cases a = 2 and 4, which corresponds to a 90° and 
180° phase lag of the xc potential with respect to the 
driving force. 

We observe that |/(2,i)| 2 strongly decreases after t w 
50 fs, after the end of the driving trapezoidal pulse. The 
opposite happens for |/(4, t)\ 2 , which starts out small and 
picks up after the driving pulse is over. Taken together, 
this is exactly what is needed to to produce the observed 
amplitude modulations of the charge-density oscillations, 
with a node approximately at t « 50 fs, see Fig. [SJ The 
amplitude suppression happens due to the contribution 
of the a = 2 TDKS determinant, where the 90° phase 
lag simulates dissipation. The subsequent increase of the 
amplitude requires the contribution of the a — 4 deter- 
minant, where the 180° phase lag simulates an elastic 
behavior. 

A similar way of explaining the amplitude modulations 
was discussed in Ref. [lj| . The time-dependent xc poten- 
tial alternatingly acts as a driving and a dissipative force, 
and is thus able to increase and reduce the amplitude of 
the charge-density oscillations. 



V. SUMMARY AND CONCLUSIONS 

The adiabatic approximation, in particular the ALDA, 
continues to be the workhorse for mainstream applica- 
tions in TDDFT, such as molecular excitation spectra. 
However, as shown by Maitra et al. [1, H|, the descrip- 
tion of more complicated processes such as multiple or 
charge-transfer excitations mandates a time-dependent 
XC functional which includes memory Since these exci- 
tation processes are extremely relevant in (bio)chemistry, 
there is an urgent need to address the question of how to 
go beyond the adiabatic approximation in TDDFT. 

In this paper, we have proposed a conceptually and 
computationally simple method to treat memory effects 
in TDDFT. Our approach is based on a new extension 
of the well-known generator coordinate (or Griffin-Hill- 
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Wheeler) method to time-dependent systems. In essence, 
we represent the A^-body wavefunction as a superposition 
of a few nonorthogonal KS Slater determinants, coming 
from a set of TDKS Hamiltonians featuring deformed XC 
potentials. The weight functions are optimized with a 
stationary action principle. The key is that our method 
allows not only spatial, but also temporal deformations 
of the XC potential. This allows for a very cheap way 
to bring in memory, namely by evaluating the XC po- 
tential at time t using a density at a previous time t r 
as input. The variational principle then determines how 
each of these retarded XC potentials are weighted rela- 
tive to each other. The numerical effort of the resulting 
GC-TDDFT scheme is quite manageable. 

We have applied our GC-TDDFT scheme to describe 
parametric oscillations in a strongly driven Hooke's atom, 
a two-electron system with the special properties that 
it cannot ionize and that its levels are nearly equally 
spaced. Therefore, the interplay between single and mul- 
tiple excitations is particularly pronounced, leading to 
a very strong modulation of the free charge-density os- 
cillations. These modulations are impossible to capture 
in ALDA or any non-retarded GC-TDDFT scheme. By 
contrast, our GC-TDDFT with memory reproduced the 
effect extremely well. 

The GC-TDDFT allows for a wide range of possible 
deformations in the TDKS Hamiltonian. As we have 
shown, one can choose deformations which remove part 
of the self-interaction error, and which bring in retarda- 
tion in a simple way. The success of the GC-TDDFT 
approach thus crucially depends on choosing deforma- 
tions according to the particular physical characteristics 
and requirements of the systems under study. We believe 
that, due to its flexibility and computational simplicity, 
the GC-TDDFT approach may prove valuable for a wide 
range of applications in chemistry and physics. 
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APPENDIX A: PROPERTIES OF THE TDGHW 
EQUATION 



The resulting static GHW equation is given by 

da' [K(a, a') - ES{a, a')} /(a') = . (A2) 



The static GHW equation provides the initial condition 
for the subsequent time propagation with the TDGHW 
scheme. From the stationary-action principle 







alt' dp d 7 /*(/3,t') 



x ($Ks(/?,0 







<i>Ks(7,0)/(7,0(A3) 



we obtain the TDGHW equation 



da 



A(a, a , t) + S(a, a , t)i — 



f(a',t) = 0, (A4) 



where the action and overlap matrix elements A and S 
are defined in eqs. (TTU)) and (jlip . 



1. Static limit 

Let us first show that the TDGHW equation reduces 
to the static GHW equation in the special case where 
the external one-body potential is time-independent. In 
that case, the KS Slater determinant associated with the 
deformation parameter a becomes 



$ks(<M) — ► e lB Ks<$ Kg ( a ) 



(A5) 



where E^ s = Ylj=i £ f ^ s tne total KS energy of the sys- 
tem. Equation (|A4[) then becomes 







da'e- iE ^ 



+ S(a,a')i^ 



E£ s S(a,a')-K(a,a') 
/(«',*), (A6) 



which reduces to the static GHW equation (|A2[) if 



/(a,t)— ►e-^-^Vfc). (A7) 



Thus, for a time-independent Hamiltonian H the weight 
functions only have a simple time-dependent phase fac- 



tor, and we obtain 



-,-iEt 



4" as expected. 



The static GHW equation follows from the varia- 
tional minimization of E = (\I>|i7|\f f )/(\I>|\I>). With the 
Hamiltonian and overlap matrix elements of the given 
DFT seed functions, K(a,a') = ($ K s(ot)\H\$Ks(a')) 
and S(a,a') — ( ( &ks(q ; )| ( I ) ks(q ! ')}; we can write this as 

8 Jdf3fd 1 r(P)K(P^)f(l) _ () fA1 s 
6f*(a) rd/3fd 7 /*(/W,7)/(7) ' ^ ' 



2. Norm conservation 

The derivation of the TDGHW equation (IA4p implic- 
itly assumes that the initial many-body wave functions 
are normalized, and that the norm 

(#(t)|*(t)>= fda[da / f*(a,t)S{a,a',t)f(a l ,t) (A8) 
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remains constant for all times. Indeed, using the defini- 
tion of S(a,a',t) and eq. (|A4[) , it is straightforward to 
show that 



dt 



(A9) 



which establishes norm conservation of the TDGHW 
scheme. 



APPENDIX B: NUMERICAL SOLUTION OF 
THE TDGWH EQUATION 

By discretizing the generator coordinate a, the GHW 
and TDGHW integral equations are turned into linear 
algebra problems. The static GHW equation becomes a 
generalized eigenvalue problem (4lj |: 



^ ^ [-^a.a' E Sa^ct'] fa 



. 



(Bl) 



where 



F = Rf. 



(B6) 



The advantage of working with the vector F is that it 
satisfies a formally simpler normalization condition com- 
pared to the original vector / of the weight functions, 
which involves the overlap matrix. It is straightforward 
to transform the TDGHW equation (|B2[) into an equa- 
tion of motion for F: 



where 



_ _ .8 

M + l di 



F = 



(B7) 



Next, we consider the TDGHW equation, where from 
now on we omit the time argument for notational sim- 
plicity: 



E 



A Q ' ® 

**-a ql' H~ &ol,ol' ^ 77~ 

d t 



fa' -0 



(B2) 



which, as we have shown above, conserves the norm 

N = J2faSaa'fa' ■ (B3) 



Calculation of the square root of the overlap matrix S, 
defined as 



(B4) 



4 R ! R 



(B8) 



idS/dt it can be shown 



M = R *AR 1 



Using Rt = R and At = A 
that Mt = M. 



The numerical solution of the equation of motion (|B7|) 
for F can now be carried out using the standard Crank- 
Nicholson algorithm [53| . which is unitary and accurate 
to second order in the time step: 



iAt 



M 



F(t + At) = 



iAt 



M 



F(t) , (B9) 



is a standard task in linear algebra. Both S and R are 
hermitian matrices. Thus, the norm (|B3[) can be rewrit- 
ten as 



(B5) 



where M is evaluated at time t + At/2. Once the time 
evolution of F has been calculated, it is a simple matter 
to obtain the time-dependent weights / = R _1 F. 
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